home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / obsolete / stepwise.pro < prev    next >
Text File  |  1997-07-08  |  12KB  |  479 lines
  1. ; $Id: stepwise.pro,v 1.2 1997/01/15 04:02:19 ali Exp $
  2. ;
  3. ;  Copyright (c) 1991-1997, Research Systems Inc.  All rights
  4. ;  reserved. Unauthorized reproduction prohibited.
  5.  
  6.  
  7.   pro writereport,X,Y,V,Enter,direction,RepType,unit
  8. ;Writereportprints out 
  9.  
  10.  printf,unit," "
  11.  printf,unit," "
  12.  printf,unit," "
  13.  printf,unit,RepType," Report"
  14.  printf,unit,direction ,Enter
  15.  printf,unit,""
  16.  regression,X,Y,VarNames=V,Unit=unit
  17.  return
  18.  end
  19.  
  20.  
  21. function setnames,vn,N,unit
  22. ; if vn is an array with n names , simply return vn,else
  23. ; pad out vn with var i.
  24.  
  25.  NC=Size(vn)
  26.  
  27.  if(NC(1) NE 0 ) THEN BEGIN        ; 
  28.    VarName=vn
  29.    if(NC(1) LT N) THEN BEGIN
  30.      printf,unit,'stepwise-missing Variable names'
  31.      I=Indgen(N)
  32.      VarName=[VarName,'Var'+STRTRIM(I(NC(1):N-1),2)]
  33.    ENDIF
  34.  ENDIF ELSE VarName='Var'+ STRTrim(INDgen(N),2)
  35.  
  36.  return,VarName
  37.  END
  38.  
  39.  
  40. function mult_cor,A,index,k
  41. return, A(k,k) - A(index,k)*A(k,index)/A(index,index)
  42. end
  43.  
  44.  
  45.  
  46. pro InsertVar,A,index,k
  47.  
  48.  mult0 = 1.0/A(index,index)
  49.  m1 = A(*,index) # replicate(1.0,k+1)
  50.  m2 = replicate(1.0,k+1)#A(index,*)
  51.  mult1 = m1*m2
  52.  case index of
  53.  0: A(1:*,1:*) = A(1:*,1:*) - mult0*mult1(1:*,1:*) 
  54.  k: A(0:k-1,0:k-1) = A(0:k-1,0:k-1) - mult0*mult1(0:k-1,0:k-1) 
  55.  ELSE:BEGIN
  56.        IN = [lindgen(index),lindgen(k-index) + index+1 ]
  57.        A(IN,0:INDEX-1) =  A(IN,0:INDEX-1) -   $
  58.                           mult0*mult1(IN,0:Index-1)
  59.        A(IN,INDEX+1:*) =  A(IN,INDEX+1:*) -  $
  60.                           mult0*mult1(IN,INDEX+1:*)
  61.       END
  62.  ENDCASE
  63.  A(index,*) = -A(index,*) * mult0
  64.  A(*,index) = A(*,index) * mult0
  65.  A(index,index) = mult0
  66. return
  67. end
  68.  
  69.  
  70.  
  71. pro UpDate, COR, INEQ, N,K
  72. ; update correlation matrix by adding variables in  $
  73. ; InEQ to equation
  74.  
  75. for i = 0L,n-1 do insertvar,COR,INEQ(i),k
  76. return
  77. end
  78.  
  79.  
  80. function remove_from, A, B,VarNo
  81. ; View both A and B as sets and return the set compliment B-A
  82. ; elements not returned in same order
  83.  
  84. Temp = fltarr(VarNo + 1)
  85. Temp(B) = 1
  86. Temp(A) = 0
  87. return, where(Temp eq 1)
  88. END  
  89.  
  90.  
  91. function remove1_from, A, B
  92. ; View both A and B as sets and return the set compliment B-A
  93. ; elements not returned in same order
  94.  
  95.  A = A(sort(A))
  96.  B = B(sort(B))
  97.  
  98.  SA = N_Elements(A)
  99.  SB = N_Elements(B)
  100.  
  101.  if SB EQ 0 THEN return,B                   $
  102.   else if SB EQ 1 THEN if A(0) EQ B(0) THEN return,C                $
  103.                        else return,B
  104.  
  105. i = 0
  106. j = 0
  107.  
  108. while i LT SA and j LT SB do BEGIN
  109.  
  110.      if A(i) EQ B(j) THEN BEGIN              ; remove B(j) from B
  111.         if j EQ 0 THEN B = B(1:*)            $
  112.          else if j EQ SB-1 THEN B = B(0:SB-2)         $
  113.               else B = [B(0:j-1), B(j+1:*)]  
  114.         i = i + 1
  115.         SB = SB - 1
  116.      ENDIF ELSE if A(i) GT B(j) THEN j = j+1  $
  117.                 else i = i+1
  118. ENDWHILE
  119.  
  120. return,B
  121. END
  122.   
  123.  
  124.  
  125.  
  126.  
  127.  
  128. function f_to_exit,C,index,COR,k,n
  129. ; f_to_exit returns the index of the variable in C with
  130. ; smallest F_to_exit value which is computed via the
  131. ; partial correlation coefficient of each variable in
  132. ;  C against Y and adjusted for the remaining
  133. ; variables in C. If C is empty, inf is returned
  134.  
  135.  SC = size(C)
  136.  SC = SC(1)
  137.  
  138.  FMin = 1.0e30
  139.  index=0
  140.  mul = n - SC - 1
  141.  
  142.  P2 = Cor(k,k)
  143.  if N_Elements(C) EQ 0 THEN return,1.0e30
  144.  
  145.  for i = 0L,SC-1 DO BEGIN
  146.     P1 = mult_cor(COR,C(i),k) 
  147.     F = (P1-P2)/P1
  148.     F = (F * mul) / (1 - F)    
  149.     if (F LT FMin) THEN BEGIN
  150.         FMin = F
  151.         index = i
  152.     ENDIF
  153. ENDFOR
  154.  
  155. return, FMIN
  156. END
  157.           
  158.  
  159.  
  160.  
  161.  function f_to_enter, X,index,COR,k,n,SC
  162. ; f_to _enter returns the maximum f_to_enter value for the
  163. ; Variables in x ,index returns the index of the var with
  164. ; max f_to_enter
  165.  
  166.  SX1 = size(X)
  167.  SX = SX1(1)
  168.  FMax = 0
  169.  mul = n-2-SC
  170.  P1 = COR(k,k)
  171.  for i = 0L, SX-1 DO BEGIN
  172.      P2 = mult_cor(Cor,X(i),k)
  173.      F = (P1-P2)/P1
  174.      if F ge 1 THEN BEGIN
  175.        index = i
  176.        return,1.0e30
  177.      ENDIF
  178.  
  179.        F = (F * mul) / (1 - F)
  180.        if (F GT FMax) THEN BEGIN
  181.           FMax = F
  182.           Index = i
  183.        ENDIF
  184.       ENDFOR
  185.  
  186. return,FMax
  187. END
  188.  
  189.   
  190.  
  191.       
  192.  pro stepwise,X1,Y1,InEQ,alpha,alphar,VarNames=VN, $
  193.                   Report=RP, List_Name =LN,Missing=M
  194. ;+
  195. ; NAME:
  196. ;    STEPWISE
  197. ;
  198. ; PURPOSE:
  199. ;    To select the locally best set of independent variables to be used
  200. ;    in a regression analysis to predict Y.
  201. ;
  202. ; CATEGORY:
  203. ;    Statistics.
  204. ;
  205. ; CALLING SEQUENCE:
  206. ;    STEPWISE, X, Y [,InEQ, Alpha, Alphar]
  207. ;
  208. ; INPUT:
  209. ;    X:    A 2-dimensional array of observed variable values. 
  210. ;        X(i,j) = the value of the ith variable in the jth observation.
  211. ;
  212. ;    Y:    A column vector of observed dependent variable values. 
  213. ;        X and Y must have the same number of rows.
  214. ;
  215. ; OPTIONAL INPUTS:
  216. ;    InEq:    Vector of indices of equations to be forced into the equation.
  217. ;
  218. ;    Alpha:    Significance level for including a variable.  The default 
  219. ;        is 0.05.
  220. ;
  221. ;    Alphar:    Significance level for removing a variable from the 
  222. ;        regression equation. The default is equal to Alpha.
  223. ;
  224. ; KEYWORDS:   
  225. ;     VARNAMES:    A vector of names for the independent variables to be used in
  226. ;        the output.
  227. ;         
  228. ;    REPORT:    A flag to print out regression statistics as each a variable 
  229. ;        enters or leaves the equation.
  230. ;
  231. ;    LIST_NAME:    A string that contains the name of output file.  Default is 
  232. ;        to the screen.
  233. ;
  234. ;      MISSING:    Value of missing data.  If this keyword is not set, assume no
  235. ;        missing data.  Listwise handling of missing data.    
  236. ;
  237. ; OUTPUT:
  238. ;    Regression statistics for the final equation.
  239. ;
  240. ; OPTIONAL OUTPUT PARAMETERS:
  241. ;    InEq:    Vector of indices of variables in final equation.
  242. ;
  243. ; RESTRICTIONS:
  244. ;    None.
  245. ;
  246. ; COMMON BLOCKS:
  247. ;    None.
  248. ;
  249. ; SIDE EFFECTS:
  250. ;    None. 
  251. ; 
  252. ; PROCEDURE:                     
  253. ;    Adapted from an algorithm from Statistical Analysis, a Computer
  254. ;    Oriented Approach, by A.A. Afifi and S.P. Azen.  The procedure 
  255. ;    successively picks the best predictor of the dependent variable Y 
  256. ;    given the variables already entered (via a statistic based on the
  257. ;    partial correlation coefficient).  Variables may be removed from the
  258. ;    equation if their significance falls below a given level. 
  259. ;-         
  260.  
  261.  On_Error,2    
  262.  
  263.  
  264.  if ( N_Elements(LN) NE 0 ) THEN openw,unit,/get,LN    $
  265.  else unit = -1 
  266.  
  267.  if N_ELEMENTS(X1) eq 0 THEN BEGIN
  268.    printf,unit,"stepwise-- Data array is empty"
  269.    goto,done0
  270.  ENDIF
  271.  
  272.  X = X1
  273.  SX = size(X)
  274.  S1 = SX(1)
  275.  S2 = SX(2)
  276.  SY = size (Y1)
  277.  
  278.  
  279.  
  280.  
  281.  
  282.  if (SY(0) LT 2) THEN         BEGIN
  283.   Y = transpose(Y1)
  284.   SY = size(Y)
  285.  ENDIF ELSE Y = Y1
  286.    
  287.  
  288.  if  SY(0) NE 2 OR SX(0) NE 2 OR S2 NE SY(2) THEN   BEGIN
  289.     printf,unit," stepwise - incompatible arrays"
  290.     goto, DONE0
  291.  ENDIF
  292.  
  293.  if N_Elements(M) THEN BEGIN
  294.    X = listwise([x,y],M,rownum,rows,here)  ; Remove cases with 
  295.                                            ; missing data.
  296.                                            ; Remove and save
  297.                                            ; corresponding
  298.                                            ; values from y
  299.    X = X(0:S1-1,*)
  300.    if N_elements(X) le 1 THEN BEGIN 
  301.       printf,unit,'stepwise---too many missing data values,'
  302.       printf,unit,'           need more observations'
  303.       Goto, DONE0
  304.    ENDIF
  305.                 
  306.    yval = Y(rownum) 
  307.    Y = Y(0,here)                          
  308.    SX = size(X)
  309.    S2 = SX(2)
  310.    SY = size (Y)
  311.  ENDIF
  312.  
  313.  COR = Correl_Matrix([X,Y])
  314.  D = Determ(COR)
  315.  
  316.  if ( abs(D) lt 1.0e-7) THEN BEGIN
  317.     printf,unit, "Determinat of correlation matrix is", D
  318.     printf,unit, " Hit control C to terminate computations" 
  319.  endif
  320.  
  321.  
  322.  
  323.  outnum = S1-1
  324.  
  325.  
  326.  
  327.  if ( N_Elements(RP) NE 0) THEN RP =1 else RP =0
  328.  
  329.  Vars = setnames(VN,S1,unit)               ;Set up variable names
  330.                                           ; for output.
  331.  
  332.  if ( N_Elements(alpha) EQ 0) THEN alpha =.05
  333.  if ( N_Elements(alphar) EQ 0) THEN alphar = alpha   
  334.  
  335.  NotinEQ = INDGEN( S1)                    ; Initially, no   
  336.                                           ; variables in    
  337.                                           ; regression eq.
  338.  
  339.  
  340.  NF = N_Elements(InEq)                  ; Get number of
  341.                                         ; variables to 
  342.                                         ; force into equation
  343.  innum = NF + 1
  344.  if ( NF EQ S1) THEN goto,done       ; If all vars are forced,
  345.                                      ; exit
  346.  
  347.  if NF NE 0  THEN BEGIN               ; Are there variables to
  348.                                       ; force  into eq$
  349.  
  350.  
  351. NotInEq = remove_from( InEq, NotInEq,S1)  ;If so,remove them 
  352.                                           ;from NotIneq list.
  353.    outnum = outnum - NF
  354. Update,COR,INEQ,NF,S1             ; Update correlation matrix
  355.  ENDIF
  356.  
  357.  
  358.  if NF NE 0 THEN C = X(InEq,*)
  359.  
  360. if (innum gt 1) THEN      $
  361.  F = f_to_enter(NotInEQ,index,COR,S1,S2,innum-1)  $
  362.                                   ;Index specifies first .
  363. ELSE BEGIN
  364.   F = max(abs(COR(S1,0:S1-1)),index)
  365.   if F ne 1 THEN F = F^2*(S2-innum-1)/(1-F^2) else F = 1.0e30
  366. END
  367.  
  368.  F = 1 - f_test1(F,1,S2-2)              
  369.                                   ;Compute significance of F.
  370.  
  371.                                          
  372.  
  373.  if ( F LT alpha) THEN BEGIN
  374.    Update,COR,NotInEq(index),1,S1
  375.    if (NF  NE 0) THEN BEGIN  ; if significant,
  376.        InEq =  [InEq,NotInEq(index)]                 
  377.    ENDIF else           $
  378.          InEq = [NotInEq(index)]
  379.  ENDIF else if (NF NE 0) THEN goto, DONE  $
  380.             else  BEGIN
  381.                  printf,unit,                                $
  382.  "stepwise- no variable significant, mean best predictor of y"
  383.       goto, DONE0;
  384.            ENDELSE 
  385.  
  386.  
  387.  if RP NE 0 THEN writereport,X(InEq,*),Y,Vars(InEq)    $
  388.             ,Vars(NotInEq(index)),     $
  389.                 "Entering Variable: ","Step", unit  
  390.  
  391.  if (outnum EQ 0) THEN goto,DONE           ;All vars in eq?
  392.  
  393.  
  394.  if ( Index EQ 0) THEN NotInEq = NotInEQ(1:*)     $
  395.                                           ; remove from list 
  396.    else if (Index EQ outnum) THEN NotInEq =     $
  397.                                NotInEq(0:outnum-1) $
  398.      ELSE NotInEq = [NotInEq(0:index-1),NotInEq(index+1:*)]
  399.  
  400.  
  401.  
  402.  
  403.  while ( innum NE S2-2 ) DO BEGIN 
  404.       F = f_to_enter(NotInEq,index,COR,S1,S2,     $
  405.                      N_Elements(InEq))      
  406.                                     ; Get next var to enter
  407.       F = 1 - f_test1(F,1,S2-2-innum)
  408.  
  409.      if ( F LT alpha) THEN BEGIN
  410.        INEQ = [Ineq,NotInEq(index)]   ; put in eq
  411.        Update,COR,NotInEQ(index),1,S1
  412.       ENDIF else goto, DONE  
  413.  
  414.      if RP NE 0 THEN writereport,X(InEq,*),Y,Vars(InEq)   $
  415.        ,Vars(NotInEq(index)),"Entering Variable: ","Step",unit
  416.  
  417.  
  418.      innum = innum +1
  419.      outnum = outnum -1 
  420.  
  421.                           
  422.      if outnum EQ 0 THEN goto,done
  423.  
  424.      if ( Index EQ 0) THEN NotInEq = NotInEQ(1:*)     $
  425.                      ; Remove it from list.
  426.  
  427.      else if (Index EQ outnum) THEN NotInEq =  $
  428.                               NotInEq(0:outnum-1) $
  429.               ELSE NotInEq = [NotInEq(0:index-1),  $
  430.                               NotInEq(index+1:*)]
  431.       
  432.      F = f_to_exit(InEq(NF:*),index,COR,S1,S2)        
  433.                         ;Are other vars still significant?
  434.  
  435.      WHILE ( 1-f_test1(F,1,S2-innum-1) GT alphar) DO BEGIN
  436.  
  437.         if NF NE 0 THEN temp = InEq(0:NF-1)  
  438.                        ; if not, remove
  439.                                         ;them from eq
  440.         innum = innum-1
  441.         outnum = outnum +1
  442.         Varout = InEq(NF + Index)
  443.         NotInEq = [NotInEQ,VarOut]
  444.         Update,COR,varout,1,S1   
  445.         if (Index EQ 0) THEN InEq = InEq(NF+1:*) $
  446.           else if (Index EQ innum -1) THEN InEq = $
  447.                                     InEq(NF:innum-2) $
  448.         else InEq = [InEq(NF:NF+index-1),InEq(NF+index+1:*)]
  449.         if NF NE 0 THEN InEQ = [temp,InEq]
  450.  
  451.         if RP NE 0 THEN writereport,X(InEq,*),Y,Vars(InEq), $
  452.           Vars( VarOut),"Exiting Variable: ","Step", unit
  453.              F = f_to_exit(InEq(NF:*),index,COR,S1,S2)
  454.      ENDWHILE
  455. ENDWHILE
  456.  
  457.      
  458. DONE:
  459.  
  460. writereport,X(InEq,*),Y,Vars(InEq),"*******","*****",  $
  461.                                      "Final",unit
  462. printf,unit," "
  463. printf,unit," "
  464. printf,unit,"Variables in Equation: ",Vars(InEq)
  465. DONE0: ;if N_Elements(M) THEN BEGIN                    ;replace cases with missing data
  466. ;   listrep,x,rownum,rows,here
  467. ;   Y = Fltarr(1,N_Elements(rownum) + sizea)
  468. ;   Y(here) = Y
  469. ;   Y(rownum) = yval
  470. ; ENDIF 
  471. if (unit NE -1) THEN Free_Lun,unit
  472. Return
  473. END
  474.  
  475.  
  476.  
  477.  
  478.  
  479.